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Abstract. We consider the Ising model for two interacting groups of spins 
embedded in an Erdos-Renyi random graph. The critical properties of the 
system are investigated by means of extensive Monte Carlo simulations. Our 
results evidence the existence of a phase transition at a value of the inter-groups 
interaction coupling which depends algebraically on the dilution of the graph 
and on the relative width of the two populations, as explained by means of 
scaling arguments. We also measure the critical exponents, which are consistent 
with those of the Curie- Weiss model, hence suggesting a wide robustness of the 
universality class. 



1. Introduction 

Statistical ferromagnetic models on diluted random graphs have received an increasing 
attention in the last few years. Although quenched randomness has been primarily 
discussed in the framework of antiferromagnetic models and competing interactions 
[U [2j [3j 0] , a renewed interest has recently emerged for quenched disorder in models 
with purely ferromagnetic interactions EJ El 13 [101 [HI 021 031 E] . Many problems 
are still open in this held, and the effects of quenched topology on the critical properties 
has not yet been exhaustively clarified. 

In d-dimensional disordered and inhomogeneous systems, the relevance of 
topological disorder at the phase transition can often be estimated by extended Harris 
criteria [15l HS1 E]> taking into account the role of geometrical fluctuations. In 
these cases, the influence of disorder on the critical behavior can be analyzed by 
estimating the shift in the critical temperature induced by disorder fluctuations in a 
region of spin variables, included in a correlated volume £ d , where £ is the correlation 
length. When the fluctuation grows in approaching the critical point, a crossover 
phenomenon can occur and the critical behavior changes. In mean field models, it is 
reasonable to expect that in high dimensions the critical behavior should not by be 
influenced by the presence of topological disorder. In this direction, one of the main 
points reached by recent studies on purely ferromagnetic models is that the mean-field 
behavior, coinciding with a complete graph topology, is very robust with respect to the 
random dilution of Erdos-Renyi random graphs [6j [18] . The universality class of the 
diluted Ising model appears independent of dilution and also topologies of graphs with 
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arbitrary degree distributions, provided the convergence of the second moment, still 
feature an analogous behavior [19 . However, a few rigorous results are available in 
this case [20, 21 , and estimates of critical exponents are often performed via numerical 
analysis. 

An interesting example of a ferromagnetic model on a random topology is 
represented by the two-populations Ising model. On the randomly diluted graph, two 
types of spins are present, with an interaction that is homogeneous among the spins 
belonging to each group, while it is characterized by a different value between the two 
groups. The possibility of considering two types of magnetic sites, interacting therefore 
through a four blocks matrix, appears an interesting point, which could be relevant 
in many phenomena, ranging from the study of anisotropic magnetic materials [22j to 
social economic models [23 , to biological systems [25]. In the latter case, an interesting 
application concerns cooperativity effects in the context of bacterial chemotaxis, where 
cooperativity between receptors has been shown to be a key mechanism for a long- 
lasting puzzle, namely how a small change in external concentration of attractant 
or repellent can cause significant amplification in receptor signal 25 , 26J . Theoretical 
models based on the two-population Ising model, although simplified, are able to show 
that a coupled system of receptors has the capacity to greatly amplify signals. These 
results suggest that the interactions among different types of receptors attain optimal 
value as it approaches the critical coupling energy of an analogous Ising model; the 
understanding of the pertaining critical properties is therefore very looked for [27] . 

Of course, all the cited applications require an extension of mean-field results 
to more general topologies. Diluted random graphs represent a first step in this 
direction. The mean field version, corresponding to a complete graph, of the two- 
populations Ising model has recently been investigated analytically [23]. There, the 
phase diagram of the model was studied, as a function of the inter-groups interactions 
J12, of the external fields and of the relative width of the two groups. In this paper we 
aim at extending these results to more general topologies, by investigating the critical 
properties of the two-population Ising model on a diluted random graph, as a function 
of the inter-groups interaction coupling J12 . First of all, we show that the simpler cases 
of two-population models on fully-connected and diluted bipartite topologies, as well 
as one-population models on Erdos-Renyi graphs can be recovered by properly tuning 
the system parameters, i.e J12, the relative width of the two groups c and the degree 
of dilution a. Moreover, non-trivial thermodynamic behavior can be achieved when 
the intra-group couplings are fixed at a small value; then, there exists a critical value 
J12 at which the system exhibits a mean-field type phase transition, in agreement 
with the fully connected model. We determine numerically the critical value of the 
inter-groups interaction coupling jf 2 , which turns out to depend algebraically on a 
and on c and we explain this algebraic behavior by means of scaling arguments. We 
also measure the critical exponents, showing that they are consistent with those of the 
Curie- Weiss model, hence suggesting a wide robustness of the universality class. 

In the following, we first describe the mathematical details of the model considered 
(Section[2j), then we present our numerical and analytical results (Section[3]). The last 
section is devoted to conclusions and final remarks (Section |4]). 
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2. The model 



Given a set of N spins denoted as {a} = {ax, a 2 , <?n}, and arranged according to 
a given topology, we define the Hamiltonian 

1 N 

H n{{o-}) = - — JijWj ~ ) ] h i a i ' (!) 



jV 

»=i 



where the first sum runs over all couples of nearest neighbours i ~ j and the whole 
set J of spins er^ is divided into two partitions 7i U I 2 = I, made up of N± = \Ii\ 
and N 2 — \I 2 \ elements, being iVj + A^2 = AT. Given two spins cr^ and Cj, their 
interaction depends on the relative subset they belong to, as encoded by the coupling 
block matrix, which reads: 



V 





N 2 






Jll 


Jl2 




} 


Ni 


J21 


J22 






> N 2 






J 







Jij > for 



In the following we consider only ferromagnetic interactions, i.e 
any i and j; also, J will be assumed to be symmetric, i.e. J 12 = J21. Therefore, 
the system under study displays three degrees of freedom: Jn and J22, tuning the 
interactions within each subgroup, and J12 tuning the interaction between spins of 
different subsets. 

As for the external field h, it assumes two values hi and /12, according to the 
subset it is applied on, as described by the following vector: 







1 


Ni 


h = 






> N 2 




\ ) 







We now introduce the magnetization for the subset Ik as 
M k ({a} k ) = -±-X2<T il 



(2) 



isl- 



and we define c = N±/N the fraction of spins belonging to Ix, it follows straightforward 
that N 2 = (1 - c)N. The total magnetization is given by M({a}) = Mx({<t}x) + 
M 2 ({a} 2 ). 

The system under investigation has been studied analytically in the special case 
of a complete graph Km (each node is linked to any other A^ — 1 nodes) as underlying 
structure [23j [24] . In this case the Hamiltonian per spin can be rewritten as 



N 



- \ [Jnc 2 M! 



hxcMx 



+ 2J 12 c(l-c)M 1 M 2 
h 2 (l-c)M 2 , 



J 22 (\-cfMl\ 



(3) 



where we dropped the obvious dependence on the magnetic configuration {a}. Then, 
assuming that states are distributed according to the Boltzmann-Gibbs measure and 
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in the thermodynamic limit, the expected value of the total magnetization was shown 
to depend on partial magnetization according to 

(AO = cMi + (1 - c)M 2 . (4) 

Moreover, the implicit analytic solution for the model reads as 

( Mi = tanh[J u cAfi + J l2 {l - c)M 2 - h x ] (r 
\ M 2 = tanh[ J X2 cM x + J 22 (l - c)Af 2 - h 2 ], ( ° } 

holding in general over the whole parameter space (c, J, h). 

Let us now turn to the case of a diluted graph as substrate. More precisely, 
we considered Erdos-Renyi random graphs Q{N,p) made up of N sites where links 
between sites are drawn independently with probability p = ot/N, in such a way that 
the average coordination number, or degree, is given by a. In this work we will always 
choose a > 1, so that the resulting graph is overpercolated. In the following we will 
neglect the external field, fixing h\ — h 2 — in order to focus on the role of a and c. 



2.1. Bounds of the model 

The Hamiltonian introduced in Eq.[7]and endowed with a symmetric J can be rewritten 
introducing the local fields (p l k , acting on the spin i and depending on the magnetic 
configuration of the subset Ik, with k = 1,2: 



Jkk 



^ = w £ ( 6 ) 



j eh 

Hence, the Hamiltonian can be rewritten as 

h n{{o-}) = -^2 ( Pi a 'i-'52 ( P2 a 'i> ( 7 ) 

ieli - g j i£li i£l2 

the first term accounts for the interaction between spins belonging to different 
populations, while the last two terms for the interaction with local fields. As long 
as the latter terms are negligible with respect to the former, the system approaches a 
bipartite, diluted Ising model. Conversely, when the latter terms prevail the system 
approaches two independent, diluted Ising ferromagnets. 

By denoting with 7 the expectation value of an observable with respect to the 
possible realizations of dilution and recalling aj = ±1, we get the following bounds 

o<KI< — — P ' ® 

where, for Nk 3> 1, the upper bounds can be approximated as Jnpc/2 and 
^22^(1 — c)/2, respectively. 

Now, whenp approaches 1, the substrate recovers the complete graph and the 
self-consistent equations ([5]) hold. Conversely, when p <C 1 (yet over the percolation 
threshold), assuming J\\/J\2 <C 1 and J22/J12 *C 1, as explained in Sec. [31 we have 
that if c is close to 1/2, the system approaches the case of a diluted ferromagnet, 
with unitary coupling, set at temperature J-^ 1 and with rescaled dilution a/2, in such 
a way that a phase transition is expected to occur at = 2/p [6]. On the other 
hand, if c <C 1, (and symmetrically if c close to 1), the smallest subset 1\ significantly 
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interacts only with spins belonging to the other group while the largest subset I2 
can display comparable interactions with both subgroups, up to very large couplings 
J12 ^ J22O-/C— 1), when we recover again the case of a diluted ferromagnet. Therefore, 
once Jn and J22 fixed as explained below, it is intuitive to look at J12 as an inverse 
"temperature" of the system, tuning the competition between an energetic and an 
entropic term: by lowering J12 a spin-flip on a node connected with both parties gets 
more and more likely. 

In the following we will focus on non-trivial cases corresponding to small c and 
small p. 

3. Numerical results 

The critical behaviour of the system described by the Hamiltonian of Eq. [7] has been 
investigated by means of Monte Carlo (MC) simulations [55] based on the single spin- 
flip Glauber algorithm according to which, given the magnetic configuration {a}, 
the spin-flip <jj — > a'j — —Oj on the j-th site, extracted randomly, is accepted with 
probability 

p({a}, a'. : , J) = . „ 77 , — r-rr , (9) 

where AHx({a}, a'j, J) = Hn({ct}, a'j) — HN{{a},aj) is the variation in the system 
energy due to spin-flip. As well-known, such a dynamic is able to lead the system 
towards equilibrium states described by the canonical distribution. More precisely, we 
checked that, once the coupling matrix J and the parameters c and a are fixed, the 
system eventually reaches a stationary state which, for N large enough so to avoid 
finite size effects, depends on the above mentioned parameters while it is independent 
of the initial configuration chosen. Moreover, the fluctuations of a specific observable 
x ({a}, a, c, J) scale like N~?. 

Therefore, the estimate of a given observable (x) is taken to be the average over 
a number of 10 3 decorrelated states of the system, once the equilibrium regime has 
been reached. The thermalization time and the decorrelation time are taken to be 
order of 10 2 MC steps; of course, as the critical region is approached (see below) the 
decorrelation time gets larger, nevertheless this choice of parameters ensures a good 
compromise between the number of measurements and the length of the run, allowing a 
proper sampling of the phase space. Moreover, the thermal averages (x) obtained this 
way are further averaged over different (~ 100) realizations of the underlying structure 
(with fixed number of nodes, average coordination number and populations width) in 
order to account for the stochasticity of the graph; however, in general, statistical 
errors due to "topological average" are significantly smaller than those arising from 
the thermal average. 

We focused in particular on the magnetizations Mi , Mi (see Eq. [2]) and on the 
susceptibility \\ the latter can be estimated as: 

X= J12KM 2 ) - (M) 2 ], (10) 

recalling that here Jn and J22 are fixed, J12 is finite and J12 S> Jn, J22 

Before proceeding we notice that we can distinguish three different cases, 
according to the values of Jn and J22: we say that Jn is small (large) if the isolated, 
i.e. J12 = 0, subsystem {ai}^^ is spontaneously in a paramagnetic (ferromagnetic) 
state; analogously for J22. Thus, if we fix Jn and J22 both large, for Jyi = we 
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have that \Mi\ = \M 2 \ = 1, and increasing J 12 the (absolute) magnetization remains 
always larger than zero. Analogously, no phase transition is expected for the case J\\ 
large (\Mi \ = 1) and J 22 small (M 2 — 0), or vice versa. In fact, by rising J 12 from 0, 
the total magnetization will increase continuously from Mi to Mi + M 2 . The effect 
is analogous to the case of a ferromagnet (system 2) set at a temperature larger than 
the critical one and coupled with an external field (system 1). 

In the following analysis we will therefore focus on the most interesting case of 
small Jn and small J 22 . 

3.1. Critical coupling J± 2 

As shown in Fig. [TJ for sufficiently large interaction strength J12 the system exhibits 
a phase transition and spontaneous magnetization occurs even in the absence of an 
external magnetic field. The value at which the phase transition occurs is denoted J^ 2 
and it is found to depend on the intra-population couplings Jn and J 22 , on the graph 
dilution a as well as on the relative size of the populations c. 




Figure 1. Main figure: Magnetization versus rescaled coupling J12 (see Eg. 1191 
here the exponent B = 1/3 was used) for systems of N = 6000 nodes, with 
p = 0.80 and c = 0.20; several choices of Jn = J22 = J are considered, as 
explained in the legend. Inset: J^ 2 as a function of J; data points (o) are fitted 
by a linear law (continuous line): in agreement with Eq. 1141 when Jn = J22 the 
critical coupling scale like ~ — J. 



In order to evaluate the critical coupling and by this the location of the transition, 
we fitted numerical data found for the magnetization and for the susceptibility using 
the following relations holding when the critical point is approached: 



M - ( 1- ^ ) . willi J vl ./;., (11) 
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and 



J 12 
Jl2 

tC 



(12) 



The measurements of Jg obtained independently from Eq. [TT] and from Eq. 
are, within the error (approximately 3%), consistent with each other. In the next 
subsection we will exploit these values to determine the critical exponents /3 and 7. 

Let us now focus on the dependence of Jg on the system parameters. As shown 
in Fig. [2j Jy2 grows as c and/or p get smaller, namely the more inhomogeneous are 
the two parties and/or the less connected is the graph, and the larger the interaction 
coupling necessary to induce a phase transition. This behavior is easy to see recalling 
that the expected number of inter-group bonds is N1N2P = Nac(l — c), so that when 
the concentration of inter-group bonds is small, the coupling strength J12 has to be 
larger in order to trigger a phase transition. 

We can further deepen this point starting from Eq.[3] As J i2 — > Jg, M\ and M 2 
are both close to zero, so that a first-order expansion leads to 



Mi w 
M 2 B 

from which we get 

{Jgf « 

namely, recalling Jn, J22 

J 12 



JncM x + Jg{l - c)M 2 , 
1 J&cMt + J 22 (l - c)M 2 



(1- J u c)[l- J 22 (l-c)] 



(1 - c)c 



«1, 
1 



(13) 



(14) 



(15) 





Figure 2. Left panel: Contour plot for measured values of J^ 2 as a function of 
p and c. Right panel: Main figure: critical coupling as a function of c for 
different values of link probability p, each represented with a different symbol, 
as shown by the legend; the continuous line represent the best fit according to 
Eq. 1151 Inset: Exponent A (see Eq. I16jl obtained by fitting data for J^ 2 versus p, 
for different values of c. Both panels refer to a system of volume N = 6000, and 
fixed internal couplings Jn = J22 = 10 — 4 . 

Let us now introduce the dilution. Our ansatz is that jg scales with p according 
to a power law, giving the following overall expression 

r P A[C) , x 

Jg « L. _ , (16) 



v^r^) 
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where the exponent A(c) depends on the inhomogeneity between the two populations; 
of course, when p = 1, the expression in Eq. IT51 found for K N is naturally restored. 

The fits performed on numerical data confirm the ansatz and show that the 
exponent A(c) approaches the value —1 from below as c goes to 0.5 (see the inset 
of Fig. [2]); otherwise stated, as the system gets more inhomogeneous, i.e. c far from 
1/2, the role of the underlying topology, encoded by p, gets more important. Moreover, 
when c is small, the dependence on c displayed by A is negligible, i.e. A(c) ~ A w —1.4 
so that the functional form of J± 2 can be factorized into J± 2 ~ P A /\/c(l — c), which 
was used to fit data in Fig. [2] (right panel, main figure). 

3.2. Critical exponents 

The power-law behaviors of Eqs. [TT] and are consistent with our numerical data 
and by means of fitting procedures we were able to estimate the critical exponents (3 
and 7, for different values of c and a. 




5 5.5 6 6.5 7 4.5 5 5.5 6 6.5 

log(|Tc-T|) log(|Tc-T|) 



Figure 3. System made up of N = 6000 nodes, with J\\ = J22 = 10 
and p = 0.80. Main figures: logM (left panel) and logx (right panel) versus 
log \Tc — T\ for the case c = 0.18. Insets: critical exponent (3 (left panel) and 7 
(right panel) for different values of c. 

In Fig. |31 we show an example of fitting curve for /3 and 7, respectively, found for 
a given realization of the substrate with c = 0.18; values of the exponents obtained in 
this way for several configurations at fixed c are then averaged in order to get the final 
estimate depicted in the inset as a function of concentrations. Our estimates of the 
critical exponents are, within the error (~ 5%), independent of c and are consistent 
with the values known for the Curie- Weiss model and the diluted Ising ferromagnet 
[6], i.e. (3 = 1/2 and 7 = — 1. Similar results were obtained also for different values of 
connectivity (a > N/4). Therefore, such a universality class exhibits a large degree 
of robustness: its properties are not affected by the (uncorrelated) dilution of the 
underlying structure and neither by a block-matrix coupling. 



A Two-populations Ising model on diluted Random Graphs 



9 



3.3. Coupling rescaling 

The expression found for J^ 2 allows to rescale the curves of the magnetization: for 
fixed J\\ = J 21 = J, we write the magnetization of each population as a function of a 
rescaled coupling J i2 

M 1 =f 1 (J 12 ,a,c), (17) 

M 2 =fc(Ji2,a,c), (18) 

where j\ and f 2 are proper functions of the system parameters and 

Jl2 = (Jl2- J?2)J B[C) - (19) 

Following the arguments and the results of Sec. 13. 1[ we expect that a large 
concentration of inter-population links (c = 1/2) gives rise to a less intensive 
dependence on J, i.e. B is close to zero, and vice versa. In fact, for c = 1/2, we 
can first consider the MF equations (Eq. [5]) and notice that, due to the symmetry, 
Mi = ta,nh(JM 1 /2 + J 12 M 1 /2) = tanh(( J 12 + 2)M 1 /2), and analogously for M 2 , where 
we used J^ 2 « 2 — J, holding for c = 1/2. Hence, Mi and M 2 both depends on J 12 
only, which yields 5(1/2) = 0. On the other hand, for the homogeneous case c = (or, 
symmetrically, c = 1), the one-population diluted Ising model with constant coupling 
J is recovered [6]. In Fig. Q]the magnetization for c = 0.20 and p = 0.80 is plotted 
versus the rescaled coupling J12 and we used _B(0.20) = 1/3. Analogous results hold 
for different concentrations c as shown in Fig. 4, where the comparison with the MF 
theory (dashed curves) is also depicted. Notice that, for both diluted and complete 
graphs, when the system is more inhomogeneous, i.e c is small, the scaling is less 
effective. 



4. Conclusions 

In this work we have analyzed the critical properties of the two-populations Ising 
model defined on Erdos-Renyi random graphs. This model features two degrees of 
inhomogeneity interplaying, namely topological dilution and a four-blocks coupling 
matrix, which make the analytical treatment definitely awkward and, on the other 
hand, make the model rather rich also in view of possible applications. In general, a 
unitary and rigorous physical description of critical phenomena in disordered systems 
still lacks and the study of further models for which there is a general agreement in 
the behavior of the corresponding pure cases is very important. In this sense, the 
two-population Ising ferromagnet is a further suitable candidate for testing the above 
predictions that has not been previously investigated in the literature. 

The system is characterized by the topological parameter a = pN, which 
represents the average number of links stemming from a node and by the number 
of nodes N making up the system, which is divided into two subsets of cardinality 
Ni and N 2 , being N 1 + N 2 = N and c = Ni/N; moreover the inter-group coupling 
is encoded by J12, while the intra-group coupling by Ju and J22, respectively. We 
focused our attention on the case of both Ju and J22 small, while J12 is tuned playing 
the role of a "temperature" . 

We ran extensive Monte Carlo simulations and we measured magnetization and 
susceptibility for the whole system evidencing the occurring of a critical point at Jf^, 
which has been shown to scale like p A ^ / yj c(l — c); the exponent A{c) approaches 
from below the value —1 as the system gets more homogeneous, namely as c — > 1/2. 
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Figure 4. Comparison between the MF solution from Eq. [3] (dashed line) and 
numerical data for the diluted case, being p = 0.8 and c = 0.1 (upper panel) and 
c = 0.5 (lower panel); different values of J are considered, as shown in the legend. 



Therefore, the role of the underlying topology gets more important the larger the 
imbalance among the two populations. 

We also measured the critical exponents /3 and 7 showing that they are consistent 
with those pertaining to the Curie- Weiss model and to the diluted Ising model. This 
result highlights the robustness of the universality class, not only with respect the 
dilution of the underlying structure, but also with respect to the existence of two 
populations. 
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